

%% in-sample period
T=50;
M=M(:,1:T);
y=y(:,1:T);
D=D(:,1:T);
X=X(:,1:T,:);


%% "true" DGP point estimates
load('../true_DGP/tDGP.mat','DGP')
true_beta=DGP(1:(2*N));


%% for Knitro

% Jacobian Sparsity Pattern:
JPattern=[repmat(eye(N),T,1),ones(N*T,2),ones(N*T,2*N),ones(N*T,N),...
    repmat(eye(N),T,1),zeros(N*T,N),ones(N*T,K),ones(N*T,KZ),eye(N*T)];
JPattern=sparse(JPattern);

% Hessian Sparsity Pattern:
% log-posterior
HPattern=zeros(6*N+2+K+KZ+N*T,6*N+2+K+KZ+N*T);
  % alpha,theta,beta,v,f
HPattern(1:(5*N+2),1:(5*N+2))=eye(5*N+2);
  % vs
HPattern(5*N+2+(1:N),:)=[zeros(N,5*N+2),eye(N),zeros(N,K+KZ),repmat(eye(N),1,T)];
  % kappa
HPattern((6*N+2+K+KZ+1):end,:)=[zeros(N*T,5*N+2),repmat(eye(N),T,1),zeros(N*T,K+KZ),eye(N*T)];
% equilibrium constraints
  % alpha
  HPattern(1:N,:)=HPattern(1:N,:)+...
        [repmat(eye(N),1,2),ones(N,2*N+2),eye(N),zeros(N),ones(N,K+KZ),repmat(eye(N),1,T)];
  % theta
  HPattern(N+(1:2),:)=HPattern(N+1,:)+ones(2,6*N+2+K+KZ+N*T);
  % f
  HPattern(4*N+(1:N),:)=HPattern(4*N+(1:N),:)+...
        [repmat(eye(N),1,2),ones(N,2*N+2),eye(N),zeros(N),ones(N,K+KZ),repmat(eye(N),1,T)];
  % beta,v,xi,gamma
for n=2:4
    HPattern((n-1)*N+2+(1:N),:)=HPattern((n-1)*N+2+(1:N),:)+...
        [ones(N,5*N+2),zeros(N),ones(N,K+KZ+N*T)];
end
HPattern(6*N+2+(1:K),:)=HPattern(6*N+2+(1:K),:)+...
        [ones(K,5*N+2),zeros(K,N),ones(K,K+KZ+N*T)];
HPattern(6*N+2+K+(1:KZ),:)=HPattern(6*N+2+K+(1:KZ),:)+...
        [ones(KZ,5*N+2),zeros(KZ,N),ones(KZ,K+KZ+N*T)];
  % kappa
HPattern((6*N+2+K+KZ+1):end,:)=HPattern((6*N+2+K+KZ+1):end,:)+...
    [repmat([repmat(eye(N),1,2),ones(N,2*N+2),eye(N),zeros(N),ones(N,K+KZ)],T,1),eye(N*T)];
HPattern=1*(HPattern>0);

HPattern=sparse(HPattern);

clear n

% Options
options=optimset('Display','off','JacobPattern',JPattern,'HessPattern',HPattern);


%% bootstrap
load('bootstrap.mat')

% % due to iterations cap imposed on Knitro (see knitro_boot.opt),
% % only "long" replication feasible here, uncomment to replicate
% NB=size(bootMAP,2);
% bootMAP=zeros(length(MAP),NB);
% exit=zeros(1,NB);
% 
% parfor bs=1:NB
%     rng(srngBoot(bs))
%     
%     % bootstrap sample
%     Dboot=zeros(N,T);
%     yboot=zeros(N,T);
%     Xboot=X;
%     % initial beliefs
%     ZZ=zeros(N,N);
%     for kz=1:KZ
%         ZZ=ZZ+Z(:,:,kz)*gamma(kz);
%     end 
%     Q=diag(v.*ep_sd)*exp(-ZZ)*diag(v.*ep_sd);
%     P=kron(eye(2),tril(Q)+tril(Q,-1)');
%     P=P\eye(2*N);
%     Pd=sqrt(diag(P\eye(2*N)));
%     B=beta;
%     for t=1:T
%         % democracy
%         kappa=mvnrnd(zeros(N,1),diag(vs.^2))';
%         U0=exp(alpha+theta(1)*(Pd(1:N,:)*sgi_nodes(:,1)'+B(1:N,:)+ep_sd*sgi_nodes(:,2)')); %#ok<*PFBNS>
%         U1=exp(alpha+theta(2)*(Pd(N+(1:N),:)*sgi_nodes(:,1)'+B(N+(1:N),:)+ep_sd*sgi_nodes(:,2)')-f-permute(X(:,t,:),[1,3,2])*xi-kappa);        
%         U0=M(:,t).*((U0./(1+U0))*sgi_weights);
%         U1=M(:,t).*((U1./(1+U1))*sgi_weights);
%         Dboot(:,t)=1*(U1>U0);
%         % growth
%         yboot(:,t)=M(:,t).*(true_beta(1:N,1).*(1-Dboot(:,t))+true_beta(N+(1:N),1).*Dboot(:,t)+mvnrnd(zeros(N,1),Sigma)');
%         % new beliefs
%         nm=find(M(:,t)==1);
%         DD=[diag(1-Dboot(nm,t)),diag(Dboot(nm,t))];
%         P([nm;N+nm],[nm;N+nm])=P([nm;N+nm],[nm;N+nm])+DD'*(Sigma(nm,nm)\DD);
%         Pd=sqrt(diag(P\eye(2*N)));
%         B([nm;N+nm])=B([nm;N+nm])+P([nm;N+nm],[nm;N+nm])\(DD'*(Sigma(nm,nm)\(yboot(nm,t)-DD*B([nm;N+nm]))));
%     end
%     
%     % boostrap estimates
%     [bootMAP(:,bs),~,exit(bs)]=knitromatlab(@(x)log_posterior(x,Dboot,M,Xboot,Z,...
%         a0,om_a,th0,om_th,b0A,b0D,om_b,s_v,d_v,f0,om_f,s_vs,d_vs),MAP,[],[],...
%         [],[],[-Inf(3*N+2,1);zeros(N,1);-Inf(N,1);zeros(N,1);-Inf(K,1);zeros(KZ,1);-Inf(N*T,1)],...
%         Inf(length(MAP),1),@(x)eq_con_jac(x,yboot,Dboot,M,Xboot,Z,sgi_nodes,sgi_weights,ep_sd,Sigma),[],...
%         options,'knitro_boot.opt');
% end


%% standard errors
se_alpha=std(bootMAP(1:N,:),0,2);
se_theta=std(bootMAP(N+(1:2),:),0,2);
se_beta=std(bootMAP(N+2+(1:2*N),:),0,2);
se_v=std(bootMAP(3*N+2+(1:N),:),0,2);
se_f=std(bootMAP(4*N+2+(1:N),:),0,2);
se_vs=std(bootMAP(5*N+2+(1:N),:),0,2);
se_xi=std(bootMAP(6*N+2+(1:K),:),0,2);
se_gamma=std(bootMAP(6*N+2+K+(1:KZ),:),0,2);


%% inference results
r_alpha=[alpha,se_alpha,2*normcdf(abs(alpha./se_alpha),'upper')];
r_theta=[theta,se_theta,2*normcdf(abs(theta./se_theta),'upper')];
r_beta=[beta,se_beta,2*normcdf(abs(beta./se_beta),'upper')];
r_v=[v,se_v,2*normcdf(abs(v./se_v),'upper')];
r_f=[f,se_f,2*normcdf(abs(f./se_f),'upper')];
r_vs=[vs,se_vs,2*normcdf(abs(vs./se_vs),'upper')];
r_xi=[xi,se_xi,2*normcdf(abs(xi./se_xi),'upper')];
r_gamma=[gamma,se_gamma,2*normcdf(abs(gamma./se_gamma),'upper')];


%% for Figure 3
results_f=table(ccodes.cnames,f,se_f);
results_f.Properties.VariableNames={'country','estimate','std_error'};
writetable(results_f,'../../Figures/Figure_3/stability_coeff_se.csv');


%% for Table A4
disp('for Table A4:')
disp(' ')
for n=1:(N-1)/2
    fprintf('%s & %5.4f & %5.4f & %5.4f & %5.4f & %5.4f & %s & %5.4f & %5.4f & %5.4f & %5.4f & %5.4f \\\\ \n',...
        ccodes.cnames{2*(n-1)+1},r_alpha(2*(n-1)+1,1),r_beta(2*(n-1)+1,1),r_beta(N+2*(n-1)+1,1),r_v(2*(n-1)+1,1),r_vs(2*(n-1)+1,1),...
        ccodes.cnames{2*(n-1)+2},r_alpha(2*(n-1)+2,1),r_beta(2*(n-1)+2,1),r_beta(N+2*(n-1)+2,1),r_v(2*(n-1)+2,1),r_vs(2*(n-1)+2,1))
    fprintf('%s & (%5.4f) & (%5.4f) & (%5.4f) & (%5.4f) & (%5.4f) & %s & (%5.4f) & (%5.4f) & (%5.4f) & (%5.4f) & (%5.4f) \\\\ \n',...
        ' ',r_alpha(2*(n-1)+1,2),r_beta(2*(n-1)+1,2),r_beta(N+2*(n-1)+1,2),r_v(2*(n-1)+1,2),r_vs(2*(n-1)+1,2),...
        ' ',r_alpha(2*(n-1)+2,2),r_beta(2*(n-1)+2,2),r_beta(N+2*(n-1)+2,2),r_v(2*(n-1)+2,2),r_vs(2*(n-1)+2,2))
end
    fprintf('%s & %5.4f & %5.4f & %5.4f & %5.4f & %5.4f &  \\\\ \n',...
        ccodes.cnames{N},r_alpha(N,1),r_beta(N,1),r_beta(N+N,1),r_v(N,1),r_vs(N,1))
    fprintf('%s & (%5.4f) & (%5.4f) & (%5.4f) & (%5.4f) & (%5.4f) &  \\\\ \n',...
        ' ',r_alpha(N,2),r_beta(N,2),r_beta(N+N,2),r_v(N,2),r_vs(N,2))
disp(' ')
disp(' ')



